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ABSTRACT 

We present the results of an X-ray mass analysis of the early-type galaxy NGC 4636, using Chandra 
data. We have compared the X-ray mass density profile with that derived from a dynamical analysis 
of the system's globular clusters (GCs). Given the observed interaction between the central active 
galactic nucleus and the X-ray emitting gas in NGC 4636, we would expect to see a discrepancy in the 
masses recovered by the two methods. Such a discrepancy exists within the central ~10kpc, which 
we interpret as the result of non-thermal pressure support or a local inflow. However, over the radial 
range ~ 10-30 kpc, the mass profiles agree within the la errors, indicating that even in this highly 
disturbed system, agreement can be sought at an acceptable level of significance over intermediate 
radii, with both methods also indicating the need for a dark matter halo. However, at radii larger than 
30 kpc, the X-ray mass exceeds the dynamical mass, by a factor of 4-5 at the largest disagreement. 
A Fully Bayesian Significance Test finds no statistical reason to reject our assumption of velocity 
isotropy, and an analysis of X-ray mass profiles in different directions from the galaxy centre suggests 
that local disturbances at large radius are not the cause of the discrepancy. We instead attribute the 
discrepancy to the paucity of GC kinematics at large radius, coupled with not knowing the overall 
state of the gas at the radius where we are reaching the group regime (>30 kpc), or a combination of 
the two. 

Subject headings: galaxies: elliptical and lenticular, cD — galaxies: individual (NGC 4636) — X-rays: 
galaxies — galaxies: kinematics and dynamics 



1. INTRODUCTION 

The current paradigm of galaxy formation describes 
how galaxies form embedded in massive dark matter ha- 
los. Whereas the measurement of rotation curves can be 
successfully applied to late-type galaxies to infer the pres- 
ence of this dark matter (see e.g. Sofue & Rubin 2001, for 
a review), this cannot be employed in early- type galax- 
ies as their stars and gas are not supported by rotation. 
Therefore, different methods must be invoked to measure 
the galaxy mass. 

It has long been known that early-type galaxies contain 
hot (^10 6 K) X-ray emitting gas (Forman et al. 1985), the 
temperature and density of which allow the determina- 
tion of the total gravitating mass, assuming hydrostatic 
equilibrium and spherical symmetry. This approach has 
proved successful at yielding meaningful mass profiles 
for early-type galaxies (e.g. O'Sullivan & Ponman 2004; 
Fukazawa et al. 2006; Humphrey et al. 2006; Zhang et al. 
2007). The effect of the assumption of spherical symme- 
try has been addressed in the case of galaxy clusters, indi- 
cating that although compression and elongation along 
the linc-of-sight can under or over-estimate the central 
mass respectively, this is only a small effect at large ra- 
dius (Piffaretti et al. 2003). However, the validity of the 
intrinsic assumption of hydrostatic equilibrium has been 
questioned with specific reference to early-type galaxies 
(Diehl & Statler 2007). NGC 4636 presents an ideal test- 
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bed in this respect, as it is a highly disturbed system, 
with evidence of bubbles and shocks caused by previ- 
ous AGN outbursts (Jones et al. 2002; Ohto et al. 2003; 
O'Sullivan et al. 2005). 

The use of dynamical tracers of the gravitational po- 
tential is also a well-established method to recover the 
kinematics of bound systems, both on the scale of glob- 
ular clusters (e.g. Gebhardt et al. 2000, in M15), and 
for the Galaxy itself (Chakrabarty & Saha 2001; Genzel 
et al. 2000; Ghez et al. 1998). The use of globular clusters 
(GCs) as tracers of the potential has been particularly 
successful in recovering mass profiles of nearby ellipti- 
cal galaxies (examples include Romanowsky & Kochanek 
2001; Cote et al. 2003; Bergond et al. 2006; Schuberth 
et al. 2006; Woodley et al. 2007). Similarly, dedicated 
surveys of planetary nebulae in early-type galaxies can 
also be used to derive the distribution of matter (Dou- 
glas et al. 2007), although care is required to avoid com- 
plications from distinct populations of planetary nebu- 
lae, which have been seen for example in the galaxy 
NGC 4697 (Sambhus et al. 2006). These approaches 
involve solving the Jeans equations under the assump- 
tion of spherical symmetry to determine the galaxy mass. 
Interestingly, a study of three early-type galaxies using 
planetary nebulae kinematics, Romanowsky et al. (2003) 
concluded a significant lack of dark matter in these sys- 
tems. However, Dekel et al. (2005) showed these data 
to be consistent with a massive dark halo when more 
radial orbits were considered. This highlights the mass- 
anisotropy degeneracy present in this approach, which 
can be broken by considering higher order velocity mo- 
ments (Lokas et al. 2007). A further systematic effect is 
the assumption of spherical symmetry. In the case where 
the galaxy is flattened along the line-of-sight, its mass 
can be under-estimated if the system is assumed to be 
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spherically symmetric (Magorrian & Ballantyne 2001). 

As both X-ray and dynamical methods have their own 
intrinsic assumptions, the most robust constraints can be 
placed on the mass profiles of early-type galaxies when 
different methods are compared. Indeed, there is cur- 
rently an emerging attempt to use different techniques in 
a complementary manner (Romanowsky et al. 2009; Chu- 
razov et al. 2008; Samurovic & Danziger 2006; Bridges 
et al. 2006); additionally, this approach improves our un- 
derstanding of the systematics involved in each method. 
Recent work by Churazov et al. (2008) explored in de- 
tail the comparison between X-ray and optically derived 
profiles for M87 and NGC 1399, finding agreement be- 
tween the methods at the 10-20% level when looking 
at the gravitational potential. However, both of these 
systems reside at the centres of clusters, M87 being the 
centre of Virgo, and NGC 1399 the centre of Fornax, and 
in both cases, the measurement of the potential is prob- 
ing the cluster potential. In so-called 'normal' elliptical 
galaxies, the situation is much less certain. For example, 
in the galaxy NGC 3379, Pellegrini & Ciotti (2006) re- 
quire an outflow of the X-ray emitting gas to bring the 
X-ray results into agreement with the optically derived 
results. Only by the study of more systems with multi- 
ple approaches will we be able to reconcile the observed 
discrepancies. This is successful on a local scale, as in- 
dividual GCs and/or planetary nebulae need to be re- 
solved, limiting the distance to which these observations 
can be made. Investigating the wider properties of the 
dark matter halos of elliptical galaxies will require tech- 
niques such as stacked lensing (e.g. Sheldon et al. 2004; 
Hoekstra et al. 2005; Klcinhcinrich et al. 2006; Koopmans 
2006; Mandelbaum et al. 2006; Ferreras et al. 2008). 

The layout of the paper is as follows. Section 2 de- 
scribes the basic properties of NGC 4636 and Section 
3 describes our method for extracting high resolution 
mass profiles from Chandra X-ray data. In Section 4 we 
present our results and comparison to the GC analysis 
of Chakrabarty & Raychaudhury (2008), and in Section 
5 we discuss the implications of our results. 

2. NGC 4636 

We have chosen to explore the properties of the galaxy 
NGC 4636 through a detailed X-ray mass analysis and 
comparison to GC data. It is a particularly interesting 
target, as the observed bubbles and shocks seen in the 
Chandra data (Jones et al. 2002; O'Sullivan et al. 2005) 
suggest departures from hydrostatic equilibrium in the 
galaxy core. The availability of dynamical data allows 
the plausibility of the assumption of hydrostatic equilib- 
rium to be explored, and the systematics of the analysis 
methods to be investigated. 

NGC 4636 is situated in a group at the edge of 
the Virgo cluster (Nolthenius 1993), of which it is the 
brightest group galaxy (Osmond & Ponman 2004). The 
group is dynamically mature, and has a virial mass of 
3.1±l.lxl0 13 M Q (Brough et al. 2006). In studies of the 
X-ray properties of the group, Rosat PSPC observations 
have shown the galaxy to have extended X-ray emis- 
sion, reaching far beyond the optical limit of the galaxy 
(Trinchicri et al. 1994). The central regions of the galaxy 
have been studied in detail using Chandra data by Jones 
et al. (2002), who identified the presence of shocks caused 
by recent AGN activity, and by O'Sullivan et al. (2005), 



TABLE 1 

Location, scale and basic properties of NGC 4636 



R.A. (J2000.0) 


12 ft 42 d 49.9 s 


Dec. (J2000.0) 


+02°41'16" 


Rcdshift 


0.003129 


Distance 


16 Mpc 


1 arcmin = 


4.7 kpc 


log L B (L s , ) 


10.47 


log L K (L KiQ ) 


11.11 


R eff (B-band)t 


1.48' 


D-25 (B-band)t 


6.03' 


log L x (ergs -1 )* 


41.59 



+ RC3 (de Vaucoulcurs et al. 1991) 
* O'Sullivan et al. (2001) 

who found evidence for gas mixing. Analysis of XMM- 
Newton data by Finoguenov et al. (2006) found features 
that include a plume of low entropy gas approximately 
10' from the centre of the system, interpreted by the au- 
thors as evidence of stripping. The X-ray mass profile of 
the system has been previously studied by Loewenstein 
& Mushotzky (2003), who found an enclosed mass at 
35 kpc of ~ 1.5x 10 12 M Q , strong evidence for a massive 
dark matter halo in this system. NGC 4636 also hosts 
a powerful AGN, with a radio power of log Li A qhz = 
21.79 (Jetha et al. 2007). 

Following Chakrabarty & Raychaudhury (2008) (here- 
after referred to as CR08) we assume a distance of 
16 Mpc for NGC 4636 throughout. The location, scale 
and basic properties of NGC 4636 are shown in Table 
1. Right ascension, declination and redshift are from the 
NASA/IPAC Extragalactic Database (NED) 4 . The K- 
band luminosity was calculated from the 2MASS K-h&nd 
magnitude (Skrutskie et al. 2006), and the i?-band mag- 
nitude was calculated from the extinction corrected Bt 
magnitude from HyperLeda (Paturel et al. 2003). The ef- 
fective radius (R e ff) is the radius enclosing half the light 
from the galaxy, and D25 is the diameter of the isophote 
describing a surface brightness of 25 mag arcsec -2 ; both 
of these parameters are quoted here for the .B-band (de 
Vaucouleurs et al. 1991). 

3. X-RAY DATA ANALYSIS 

Our aim in this analysis is to produce a high resolution 
mass profile from Chandra data. We have achieved this 
through a two stage approach, by firstly concentrating 
on constraining the temperature profile, followed by de- 
termining the gas density profile in much greater detail. 
This is a similar technique to that employed by Vikhlinin 
et al. (2006), although here we use the xspec projct 
model to deproject the spectra. The full analysis proce- 
dure is described in detail below, and this procedure has 
already been applied to the early-type galaxy NGC 1407 
(Romanowsky et al. 2009). 

3.1. Initial data reduction 

NGC 4636 was observed by Chandra on February 14th 
2003 (obs ID = 3926) for 75.69 ks, and we used this 
archival data in the following. The initial data reduc- 
tion was performed using version 3.4 of the Chandra In- 
teractive Analysis of Observations 5 (ciao) with CALDB 

4 http://nedwww.ipac.caltech.edu/ 

5 http://cxc.harvard.edu/ciao 
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Fig. 1. — A smoothed Chandra image of NGC 4636, extracted 
across the energy range 0.3-2.0 keV. The circles show radii of 100", 
200", 300", 400" and 500" for reference, centred on the galaxy co- 
ordinates from NED. 

version 3.4.2. We extracted a new level 2 events file 
from the level 1 events file, and removed events with 
ASCA grades of 1, 5 and 7. Bad pixels were also re- 
moved from the analysis, and the appropriate gain file 
and time-dependent gain correction were applied. Flares 
were eliminated from the events file by extracting a light 
curve from each of the back-illuminated chips, and one 
from the front-illuminated chips. The light curves were 
filtered using the 'lc_clean' script of Markevitch 6 , result- 
ing in a cleaned exposure time of 74.7 ks. 

Point sources were detected using the CIAO tool 
'wavdetect', and were excluded from further analysis. 
Spectra and response files were extracted from the 
cleaned events file, following the CIAO analysis threads. 
The blank sky background files 7 were used to extract 
background spectra; this was done over the same area as 
the source regions. The background spectra were nor- 
malised at high energies (9-12 keV) to match the source 
spectra. The use of blank sky backgrounds has been 
questioned in analyses of diffuse emission (e.g. Humphrey 
et al. 2006) , but this is particularly a problem in systems 
or regions of low surface brightness, where decompos- 
ing the diffuse emission and the background components 
can be very difficult. In the case of NGC 4636, the ob- 
servation is heavily source dominated, so this is not a 
significant issue. However, we investigate the implica- 
tions of using the blank sky backgrounds on our analysis 
in Section 4.1.3. 

To indicate the scale of the disturbances in the centre 
of NGC 4636, Figure 1 shows an image extracted from 
the Chandra data, across the energy range 0.3-2.0 keV, 
and smoothed using the CIAO tool 'aconvolve'. 

3.2. Spectral analysis 

Our two-stage spectral analysis is designed to extract 
high resolution mass profiles. We initially extract spectra 
from a series of wide concentric annuli, which we term the 

6 http://cxc. harvard. edu/ciao3.4/ahclp/lc_clean. html 

7 http://cxc.harvard.edu/contrib/maxim/acisbg/ 



coarse stage, followed by extracting spectra from much 
thinner annuli during the fine stage. The coarse stage ro- 
bustly constrains the temperature, and the fine stage in- 
corporates these constraints in determining the gas den- 
sity. The procedure is explained in detail below. We 
performed all the spectral analysis using xspec Version 
11.3.2t, and all spectra were fitted in the energy range 
0.7-7.0 keV. 

3.2.1. Coarse spectral analysis 

Initially spectra and their associated responses were 
extracted from annuli centred on the galaxy co-ordinates 
from NED (shown in Table 1), and background spectra 
were extracted from the blank sky backgrounds. The an- 
nuli were chosen to contain a net number of counts that 
allowed for both a successful deprojection of the spectra, 
as well as placing robust constraints on the temperature. 
It was found that a criterion of 8000 net counts per spec- 
trum was more than adequate, yielding 13 radial bins. 
This criterion could be further relaxed and still provide 
a good fit to the spectra, however due to the nature of our 
approach, we probe the gas on a finer radial scale in the 
second stage of the analysis. The central 0.3' (~1.4 kpc) 
was excluded from the analysis, due to the sudden peak 
in surface brightness in the image in this region. 

We fitted absorbed APEC models in each annulus, 
with an additional power-law component subject to the 
same absorption to constrain the contribution from un- 
resolved Low Mass X-ray Binaries (LMXBs). The ad- 
dition of this model component is explained in detail in 
Section 3.3, but it is important to note here that this 
component is modelled as a background component, and 
is not deprojected. We fixed the absorption (N-r) at 
the Galactic value of 1.82xl0 20 cm~ 2 (Dickey & Lock- 
man 1990) throughout, and all abundances are quoted 
as those of Grevesse & Sauval (1998). The spectra were 
then deprojected using the PROJCT model in xspec, un- 
der the assumption of spherical symmetry. During this 
procedure, the abundance was tied between all annuli, 
as otherwise it was unconstrained in some spectra. We 
discuss this assumption in full in Section 4.1.1. We de- 
fine a characteristic radius, r, for each annulus using the 
emission-weighted calculation of McLaughlin (1999), 

^[0.5(4 /2 +^ 2 )] 2 /3 (1) 

where r in and r out are the inner and outer radial bounds 
of the annulus respectively. The deprojection therefore 
yields three-dimensional temperature, abundance and 
Nh as a function of radius. 

We fit smoothing spline functions to the deprojected 
profiles, using the SMOOTH. SPLINE function from the R 
PROJECT statistical package (R Development Core Team 
2008), to give a functional form for each deprojected pro- 
file. In this case the abundance and Ah are constant 
as a function of radius. The benefits of the smoothing 
spline function are that it responds to natural variation in 
the profile, without imposing a prescribed analytic form. 
Statistical fluctuations are limited by weighting the fit 
using the inverse variances of the parameters from a se- 
ries of Monte Carlo realisations, which we describe in 
Section 3.5. The end product at this stage is a smooth, 
continuous functional form for the deprojected temper- 
ature profile; the data and associated smoothing spline 
function are shown in Figure 2. 
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Fig. 2. — The deprojected temperature profile for NGC 4636 and 
the associated smoothing spline fit, shown as a solid line (see text 
for details). The vertical error bars are Icr errors from 200 Monte 
Carlo realisations of the procedure and horizontal error bars show 
the radial extent of each annulus. 

3.3. Low mass X-ray binary component 

In addition to the emission from the diffuse hot gas, 
there is a hard X-ray contribution to the spectrum from 
unresolved Low Mass X-ray Binaries (LMXBs) . This was 
modelled using a power-law model component subject to 
the same absorption as the APEC model, the index of 
which was fixed at 1.56 (Irwin ct al. 2003). Although 
~50% of the LMXBs in NGC 4636 are associated with 
GCs, and the light profile of GCs differs compared to the 
halo light profile (Kim et al. 2006), the distribution of 
LMXBs is comparable to the halo light profile (Kim et al. 
2006). In treating the LMXB component, we therefore 
make the assumption that the distribution of LMXBs fol- 
lows the halo light of the galaxy, which we approximate 
with a de Vaucouleurs profile with R e ff = 1.48' (de Vau- 
couleurs et al. 1991), equal to 6.9 kpc for our assumed 
distance of 16 Mpc. The de Vaucouleurs profile is an 
appropriate choice for galaxies with R e ff greater than 
6.3 kpc (Prugniel & Simien 1997). 

We fit absorbed apec+powerlaw models to the spectra 
in xspec, to determine the normalisation of the power- 
law model component, under the constraint that the nor- 
malisation in each annulus should follow the overall shape 
of the light profile. In practice, this means that the model 
normalisations are tied in proportion to the shape of the 
light profile. We use the XSPEC command fakeit to fake a 
spectrum corresponding to the power-law model in each 
annulus, which is then added to the background spec- 
trum. This means that this model component is not 
deprojected. The same approach is applied in the fine 
stage to quantify the LMXB component. 

The contribution from LMXBs depends on the op- 
tical luminosity of the galaxy (O'Sullivan et al. 2001; 
Kim & Fabbiano 2004), allowing a consistency check 
on our adopted approach. Over the radial range cov- 
ered by our coarse bins, and in the 0.5-0.7 keV energy 



range, the total flux from the two-dimensional spectral 
fitting to the coarse spectra is 9.14xl0 -12 erg s _1 cm~ 2 . 
The flux from just the power-law model component 
is 1.49xl0~ 13 ergs^ 1 cm -2 , giving a fractional contri- 
bution to the total flux from unresolved LMXBs of 
1.6%. Assuming that the X-ray luminosity from discrete 
sources is log Ld. SC r = 29.5 erg s _1 Lg^ (O'Sullivan et al. 
2001), we expect, given a _B-band luminosity of 10.47 
(Table 1), a contribution from the discrete sources of ap- 
proximately log Ldscr — 39.97 erg s _1 . Assuming log L x 
= 41.59 O'Sullivan et al. (2001), the expected unresolved 
source contribution to the total luminosity is ~2%. Our 
unresolved flux is reasonable given this prediction. The 
high resolution of Chandra and the deep observation of 
NGC 4636 will have allowed more of the brightest point 
sources to be detected and excluded than in the Rosat 
data of O'Sullivan et al. (2001), which could easily lead 
to the slight difference in the predicted and observed un- 
resolved source fraction. 

As a thermal bremsstrahlung component with a fixed 
temperature of 7.3 keV can also be used to describe the 
LMXB spectrum (Irwin et al. 2003), we tested the use of 
this component instead of the power-law described above. 
There was no improvement in the fit, the fitted parame- 
ters were consistent with those recovered from using the 
power-law, and the percentage of the total flux in the 
0.5-7.0 keV energy range was found to be 1.5%, again 
consistent with that recovered from the fitting using the 
power-law model. 

3.4. Fine spectral analysis 

We next determine a set of finely spaced annuli, from 
which source spectra, background spectra, and the ap- 
propriate response files are extracted. Our motivation 
here is to model only the gas density; the remaining pa- 
rameters in our model are described by the fitted func- 
tions from the coarse stage. The fine annuli are spaced 
using a net counts criterion, however, as we are now only 
fitting for one parameter, the number of net counts in 
each annulus can be considerably reduced. We use 2000 
net counts per spectrum (51 spectra) as a compromise 
between resolution and the time taken to perform the 
Monte Carlo error analysis. Under this criterion, the an- 
nular width of the bins ranges between a minimum of 
2.95" and a maximum of 66.9", meaning that the bin 
width exceeds the PSF at all radii. We determined the 
characteristic fine radii for these annuli using Equation 1 . 

Using the functional fits to the deprojected profiles de- 
scribed in Section 3.2.1, we interpolated the values of de- 
projected temperature, abundance and -/Vh at the char- 
acteristic fine radii. These parameters were kept fixed 
in the subsequent PROJCT model fit to the fine spectra. 
The contribution from LMXBs was included as described 
in Section 3.3, and to speed up the fitting in this stage, 
we only extracted spectra across 8 channels. There is 
therefore just one free parameter at this stage in the de- 
projection — the APEC model normalisation, K, from 
which the gas density can be directly determined as 

10~ 14 f 

K = MD A (l + zW J nHUedV ' (2) 

where Da is the angular diameter distance to the galaxy, 
z is the redshift of the galaxy, and ur and n e are the 
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Fig. 3. — The gas density profile for NGC 4636, determined from 
the finely binned spectral analysis. The solid line is a /3-modcl fit 
to the data points, and the vertical error bars show the Icr errors 
from 200 Monte Carlo realisations of the procedure. Horizontal 
error bars show the radial extent of each annulus. The dashed 
line shows the effect on the gas density profile of allowing for the 
observed projected abundance gradient (see Section 4.1.1), and the 
associated 2cr confidence region is derived from /3-model fits to 100 
Monte Carlo realisations. 



number density of hydrogen and electrons respectively, 
and we assume that nn/n e = 1.17. Therefore, for a 
particular spherical shell of volume dV, the gas density 
can be recovered. The resulting gas density profile is 
shown in Figure 3, and has been fitted with a /3-model, 
yielding parameters of r core = 31.6" and (3 — 0.5 for a 
reduced x 2 of 4.7 (48 degrees of freedom). The benefit 
of this two stage approach is that we retain a robust 
temperature profile, but improve the resolution of the 
gas density profile. 

It is apparent from Figure 3 that the /3-model shape is 
not successful in describing the shape of the gas density 
profile at all radii, and fluctuations from this smooth pro- 
file can be seen. However, the calculation of the resulting 
mass profile requires a smooth and continuous function, 
and we shall proceed with the use of the fitted /3-model 
in this context. The observed fluctuations can be readily 
understood in terms of the observed disturbances in the 
X-ray emitting gas. Examining the X-ray image shows 
that the edge of the central shock region occurs at a 
radius of approximately 100", which corresponds to a 
slight depression in the gas density profile shown in Fig- 
ure 3. The feature at ~60" arises as a consequence of 
the shocks in the galaxy core. The nature of the spectral 
analysis and deprojection means that the profiles shown 
here represent azimuthally averaged measurements, so 
very localised features in the hot gas would be smoothed 
out. In the context of the gas density behaviour beyond 
^500", we note that Trinchieri et al. (1994) reported a 
flattening of the gas density distribution at radii of 6'-8' 
on the basis of Rosat data. We will examine this feature 
in further detail in Section 5.3.2. 



3.5. Monte Carlo error analysis 

We have employed a Monte Carlo (MC) approach to 
calculate the errors associated with the procedure, im- 
plemented in both the coarse and fine stages in an anal- 
ogous way. Initially, the best-fitting PROJCT model from 
the coarse analysis is used to produce a series of spec- 
tra using the xspec command fakeit with the inclusion 
of random Poisson noise. These spectra are then fitted 
with a PROJCT model, and the process is repeated 200 
times, from which the standard deviation is used to de- 
fine la errors on the coarse profiles. The errors are also 
calculated for the fine stage of the analysis by using the 
best-fitting PROJCT model from the fine analysis to fake 
a series of spectra. These are fitted with PROJCT models 
to determine the APEC model normalisation, from which 
la errors are determined at the fine radii. The MC re- 
alisations of the temperature profile are fitted with the 
SMOOTH. spline algorithm to determine continuous func- 
tions, which are used in conjunction with /3-model fits to 
the MC gas density profiles to yield 200 MC realisations 
of the mass profile. From this suite of mass profiles, er- 
rors are estimated by determining the la spread in the 
functions evaluated at the characteristic radii. 

4. RESULTS 

Here we present our X-ray derived mass profile, before 
comparing this to the results derived from a dynamical 
analysis of the GC population, performed by CR08. 

4.1. X-ray mass profile 

The resulting fits to the temperature (Section 3.2.1) 
and gas density profiles (Section 3.4) are used to deter- 
mine the mass within a given radius M(< r), in the 
following way (Fabricant et al. 1980), 

where p is the gas density, T is the temperature, G is the 
gravitational constant, /x is the mean molecular mass (as- 
sumed here to be 0.593 for a fully ionised plasma) and 
m p is the mass of a proton. The mass profile derived 
from the X-ray analysis is shown in Figure 4. The con- 
fidence region shows the 2a spread from the 200 Monte 
Carlo realisations of the mass profile. The X-ray gas den- 
sity and temperature measurements within 30" are well- 
constrained in the spectral fits, and are not obviously bi- 
ased. However, the calculation of the mass profile yields 
an unphysical negative mass in this region, indicating the 
requirement for additional non-thermal pressure support 
within the central ~ 3 kpc. 

Performing a mass analysis on an earlier Chandra 
dataset, Loewenstein & Mushotzky (2003) determined a 
total mass of ^1.5xl0 12 M Q at ^35 kpc. Correcting for 
the different assumed distance, we plot the enclosed mass 
recovered by Loewenstein & Mushotzky (2003) in Figure 
4. This falls just below the 2 a confidence bound on our 
original mass profile. Loewenstein & Mushotzky (2003) 
also found the mass to increase as r 12 over the radial 
range studied (0.7-35 kpc). Fitting a powerlaw model 
to our data, we find the same slope (1.19±0.01) outside 
5 kpc. Within this radius, our mass profile falls away 
more steeply, a consequence of the steeper temperature 
profile recovered from our analysis in the inner regions. 
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Fig. 4.— The cumulative mass profile of NGC 4636. The white 
dotted line and the associated black (2 a) confidence region show 
the results of our X-ray mass analysis. The confidence region has 
been determined from 200 Monte Carlo realisations of our proce- 
dure. The dashed line shows the effect on the calculated mass of 
allowing for the abundance gradient (see Section 4.1.1), and the as- 
sociated lcr confidence region shows the results of 100 Monte Carlo 
realisations of the procedure. The black dotted line shows the 
powerlaw slope of r 1,2 determined by Loewenstein & Mushotzky 
(2003), with arbitrary normalisation, and the cross point shows 
the total mass measured by Loewenstein & Mushotzky (2003) at 
35 kpc (corrected for our assumed distance). 

Figure 5 shows the total mass density profile from the 
X-ray procedure evaluated at the finely spaced radii, and 
the total mass density profile recovered when the abun- 
dance gradient is included in the fitting as detailed in 
Section 4.1.1. The mass density profile of elliptical galax- 
ies is the combination of the stellar mass density and the 
underlying dark matter density, and the stellar mass den- 
sity dominates within approximately li? e // (Mamon & 
Lokas 2005). 

4.1.1. Implications of allowing for the abundance gradient 

Recent work by Rasmussen & Ponman (2007) has 
shown abundance gradients to be prevalent in galaxy 
groups, and we now consider the implications of implic- 
itly assuming a flat abundance profile, when the pro- 
jected analysis reveals an abundance gradient (Figure 
6). In the innermost bins, the abundance is poorly con- 
strained and reaches the default xspec fitting limits - 
this was our main motivation for imposing a flat profile. 
This is probably a consequence of multiple temperature 
and abundance components in the inner regions (see Sec- 
tion 4.1.2). The fitted deprojected abundance is also 
shown in Figure 6, and it can clearly be seen that impos- 
ing this criterion underestimates the abundance within 
200" and overestimates the abundance outside this ra- 
dius. We are motivated to test the effect on the mass 
profile of assuming a constant abundance as this is often 
employed to satisfactorily constrain model parameters in 
less luminous systems, or in cases of poorer data qual- 
ity, and the effect of such an assumption has not been 
studied. It is a particularly important issue in low tem- 
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Fig. 5.— The total mass density profile of NGC 4636 derived 
from the cumulative mass profile shown in Figure 4, and evaluated 
at the characteristic fine radii. The horizontal error bars show the 
radial width of each bin, and the vertical error bars show the lcr 
spread of 200 Monte Carlo realisations evaluated at each radius. 
The dashed line shows the total mass density profile when an al- 
lowance is made for the observed abundance gradient (see Section 
4.1.1 for details). 

perature systems where line emission dominates. 

We repeated our coarse deprojection, fixing the depro- 
jected abundances at their projected values. It is very 
likely therefore, that we have now overestimated the very 
central abundance due to its poor constraints, so we con- 
sider the following to be an upper limit on the effects of 
allowing for the abundance gradient. The abundance and 
APEC model normalisation play off against each other 
due to line emission dominating the flux at low temper- 
atures, so to see the full effects on the gas density we 
proceeded with the fine stage. We fixed the temperature 
profile in the fine stage deprojection at the values interpo- 
lated from the original fit to the data (see Section 3.2.1), 
as the variation in temperature caused by allowing for 
the abundance gradient was well within the 1 a errors of 
the original temperature profile. To establish the errors 
in this analysis, we performed 100 Monte Carlo realisa- 
tions of the procedure, using the MC realisations from 
the coarse stage to weight the smoothing spline fit to the 
abundance profile at the beginning of the fine stage. 

The effect of allowing for the abundance gradient is to 
flatten the gas density at all radii. Fitting a /3-model, 
weighted by the inverse variance from the MC realisa- 
tions, gives [3 = 0.3 and r core = 27.2" and is shown in 
Figure 3. The associated lcr confidence region shows 
the range of /3-model fits allowed by the MC realisations. 
The subsequent effects on the mass profile and total den- 
sity profile are shown in Figures 4 and 5. The effect of 
the abundance gradient is to reduce the mass at all radii 
by a factor of ~1.6. This demonstrates the intricacies 
involved in the detailed X-ray analysis, as this effect is 
not allowed for by our Monte Carlo procedure, which also 
ties the abundances in the deprojected fit and hence this 
is a key systematic in the application of our method. 
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Fig. 6. — The projected abundance profile of NGC 4636 (data 
points), with errors estimated from XSPEC. In the two innermost ra- 
dial bins, the abundance is unconstrained, shown by the data points 
reaching the XSPEC default fitting limit of 5.0Zq. The dashed line 
shows the fitted abundance in the deprojection, where the abun- 
dance is tied between all the radial bins. 

4.1.2. Implications of assumed single temperature model 

Throughout our analysis we have implicitly assumed 
that the gas in each annulus is single phase. If the gas is 
multi-phase in these regions, the recovered abundances 
may be affected by "Fe-bias" (Buote 2000), where the 
abundance of multi-phase gas is underestimated if a sin- 
gle temperature model is fitted. The disturbed nature 
of the gas in NGC 4636 (Jones et al. 2002; O'Sullivan 
et al. 2005) suggests multiple temperatures and abun- 
dance will be present in each coarse region. Such an in- 
tegrated spectrum would have a broader iron peak, and 
would require very high quality spectra to separate the 
individual components. Using the 2-d coarse spectra, 
we tested the addition of an extra APEC model compo- 
nent, but found that this did not improve the fit at any 
radius. Although this system has been shown to host 
cavities (Ohto et al. 2003) and also shows surface bright- 
ness features (O'Sullivan et al. 2005), it appears that 
when considering an azimuthally averaged profile with a 
sufficiently large number of counts, a single temperature 
model is acceptable. In terms of our mass analysis, it is 
important to have a good representation of the temper- 
ature profile, even if the model itself does not give the 
most statistically accurate fit. 

4.1.3. Implications of using blank sky backgrounds 

To test the sensitivity of our results to the use of the 
blank sky backgrounds, we performed the following tests. 
Using the outermost coarse annulus, which will be the 
most sensitive to the background, we fitted a simple ab- 
sorbed APEC model with ]Vh fixed at the Galactic value, 
to recover the temperature, abundance and APEC model 
normalisation shown in Table 2. We scaled the exposure 
time of the blank sky background spectrum for this an- 
nulus up and down by 10 percent, re- fitting each time 



TABLE 2 

The recovered parameters from tests carried out on the 
use of the blank sky backgrounds. 



(keY) 



Abundance 

(go) 



APEC norm 



Blank sky 
Blank sky + 10% 
Blank sky - 10% 



0.79±0.01 
0.79±0.01 
0.79±0.01 



0.17 1 



7+0.03 

-0.02 

0.15±0.02 



1.43 



0.21 



+0.04 
-0.03 



,+0.12 
"-0.11 
+0.12 
-0.11 
,+0.11 



xlO 



Notes: The fitted model parameters are shown for an absorbed 
APEC model fit to the outermost coarse annulus, for the assumed 
blank sky background, and also for variations of ±10% in the nor- 
malisation of this background. Errors are derived from XSPEC and 
arc la. 

to effectively alter the background normalisation, as the 
same number of counts are collected over a differing time 
period. The recovered parameters are also shown in Ta- 
ble 2. We find that the fitted parameters for the increased 
exposure time are consistent within 1 a, as is the recov- 
ered temperature for the decreased exposure time case, 
with the abundance and model normalisation consistent 
with the original fit within 2 a. These tests indicate that 
our results are not sensitive to variations in the back- 
ground level at the level of 10 percent, and emphasises 
how robust the temperature measurement is at these low 
temperatures due to the dominance of the line emission. 

4.1.4. Implications of assumed LMXB model 

The i-T-band is a better description of the older stellar 
populations of early-type galaxies than the -B-band, and 
therefore may better describe the LMXB distribution. 
We considered the effects of assuming a de Vaucouleurs 
profile for the LMXB population with a if-band R e ff of 
56.2" (Jarrett et al. 2003). This assumption reduces the 
gas density at all radii by approximately 3 %, and makes 
no discernible change to the innermost nine temperature 
points (less than 1%). Instability from the deprojection 
procedure is however visible at the largest radii. The 
fitted abundance, which was again tied between the an- 
nuli, was 0.85Z Q using the if-band R e ff, compared to 
0.79Z Q in the original analysis. We therefore assert that 
the effect on the mass profile of the treatment of any 
abundance gradient is more crucial in this case than the 
intricacies of the treatment of the LMXBs. This may 
not be the case for galaxies where the unresolved source 
emission is a higher fraction of the overall X-ray emission. 

4.2. Comparison to dynamical mass estimate 

A sample of 174 GCs in NGC 4636 were tracked for 
their line-of-sight velocities in the observational pro- 
gramme of Dirsch et al. (2005) and were used to assess 
the mass profile of NGC 4636 by Schuberth et al. (2006). 
CR08 input these kinematic data into the Bayesian non- 
parametric algorithm CHASSIS (Chakrabarty & Saha 
2001). This invokes a Markov Chain Monte Carlo 
(MCMC) optimiser to recover the most likely equilib- 
rium distribution function from which this kinematic 
data could have been drawn, given the recovered poten- 
tial in which the sample of GCs resides. This potential is 
expected to be the gravitational potential of the galaxy 
itself, from which the total (luminous+dark) matter den- 
sity of NGC 4636 is estimated. Motivated by a desire to 
understand and test the underlying assumptions of these 
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Fig. 7.— The total mass density profile of NGC 4636. Solid 
squares show the results from the X-ray analysis presented here, 
with vertical error bars showing la errors from 200 Monte Carlo 
realisations of the procedure. Stars show the results from the GC 
analysis of CR08, where vertical error bars are the la spread in 
mass models derived from the CHASSIS algorithm. Note that the 
radial range of the X-ray measurements has been restricted to that 
determined by the dynamical mass measurements, and these radii 
have been used to evaluate the X-ray profile. The dashed line shows 
the recovered mass density profile when the abundances are fixed 
at their projected values (see Section 4.1.1). The dotted line and 
confidence region show the estimated stellar mass density (see text 
for details). 

two independent methods, we can view the X-ray mass 
profile in comparison to the dynamical estimate of the 
total mass distribution. At this point, it merits mention 
that in its current form, CHASSIS assumes isotropy in 
phase space, although work is underway to relax the re- 
quirement of velocity isotropy (Chakrabarty & Saha, in 
preparation) . 

Figure 7 shows the total mass density profile of 
NGC 4636 (star symbols) recovered in RUN I of CHAS- 
SIS by CR08; we refer the reader to this work for more 
information. The mass density from the X-ray analy- 
sis is shown as solid squares, and has been evaluated at 
the radii of the dynamical mass profile. The errors on 
the dynamical mass estimate indicate the ±lcr spread 
in the mass models about the most likely configuration, 
as determined by the MCMC optimiser that is used in 
CHASSIS. We show the stellar mass density determined 
from the iC-band luminosity density profile presented by 
CR08, assuming a stellar mass-to-light ratio in the K- 
band of 0.83. This is a colour-dependent estimate using 
the total B — V colour from HyperLeda 8 of 0.94, and 
converting to the if -band mass-to-light ratio following 
the prescription of Bell ct al. (2003). 

The general nature of the dynamical and X-ray mass 
distributions is similar to about 30 kpc, beyond which, 
the X-ray mass exceeds the dynamical mass, by a fac- 
tor of ^4.5 at 40 kpc. There is also an indication of a 
'break' in the dynamical mass density profile, at a little 

8 http://leda.univ-lyonl.fr/ 



over 30 kpc. Of the 174 GCs studied by Schubcrth et al. 

(2006) , only 15 of these are at radii greater than 7.5'. 
In terms of the GC density distribution, a steepening is 
observed between approximately 6' and 8' (Dirsch et al. 
2005), noted by Schuberth et al. (2006) to be inconsistent 
with NGC 4636 being in a dark matter potential which 
smoothly reaches to large radius. We will return to this 
point in Section 5. The key difference between the pro- 
files at large radii is the shape; the GC profile appears 
to 'break' at approximately 400", dropping away more 
steeply than the X-ray derived profiles. We note that 
the effect of allowing for the observed abundance gradi- 
ent in the X-ray analysis reduces the mass density at all 
radii, improving the agreement in the outer regions. 

CR08 fit a Navarro et al. (1996) density profile, 
hereafter NFW profile, to the total density outside 
~32 kpc recovered from the GC analysis using CHAS- 
SIS. They find a concentration of 9, and scale radius r s of 
33.7 kpc±ll percent. We fit the X-ray mass density pro- 
file across the radial range shown in Figure 7 using the R 
PROJECT non-linear least squares algorithm 'nls', weight- 
ing each point by its inverse variance. We find a concen- 
tration of 20.1±0.8 and a scale radius of 21.8±0.9 kpc, 
where the quoted errors are 1 a standard errors on the 
fit. However, this concentration is an over-estimate due 
to ignoring the stellar contribution to the mass density, 
so we proceed to fit an NFW profile to the mass density, 
having subtracted the stellar mass density shown in Fig- 
ure 7. There is some uncertainty in this approach due to 
our assumed stellar mass-to-light ratio, but this fit does 
recover a lower concentration (18.0±0.6), with a scale 
radius of 24.6±0.9 kpc, leading to an estimate for r2oo 
of approximately 443 kpc. This concentration is similar 
to NGC 720 and NGC 1407 (Buote et al. 2007), which 
are slightly cooler and warmer (~0.5 keV and ~1.0 keV; 
Osmond & Ponman 2004) than NGC 4636 respectively 
Further increasing the stellar component by increasing 
the mass-to-light ratio would further reduce the recov- 
ered concentration. 

If we instead fit the NFW profile to the X-ray profile 
where the abundance gradient has been incorporated, we 
find a concentration of 14.4±0.4 when we first subtract 
the stellar mass, with a scale radius of 26.1±1 kpc re- 
spectively. We show the NFW fits to the X-ray analysis 
when the stellar mass is subtracted in Figure 8. These 
solutions also both fit on the c-M relation of Buote et al. 

(2007) . The shape of the NFW profile is a good approx- 
imation to the shape of the recovered total mass density 
profile, and it would be difficult to fit an NFW profile 
across the whole radial range to the GC data, due to the 
small 'break' in the profile at radii of ~400-500". 



4.2.1. Comparison of gravitational potential 

Perhaps a more appropriate method of comparison is 
to look at the gravitational potential recovered from each 
method. As shown by Churazov et al. (2008), the grav- 
itational potential can be easily recovered from an X- 
ray analysis of this type, if hydrostatic equilibrium and 
spherical symmetry are further assumed. The require- 
ment that the thermal gas pressure is the only contribu- 
tor to the overall pressure, yields the following expression 
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Fig. 8.— The total mass density profile for NGC 4636. The data 
points are as shown in Figure 7. The solid line is the NFW fit to the 
main X-ray analysis, having subtracted the stellar mass component 
shown in Figure 7. The dotted line is the NFW fit to the X-ray 
analysis when the abundance gradient has been allowed for, and 
again the stellar mass component was subtracted to perform the 
fit. The short, vertical solid and dotted lines on the x axis are the 
scale radii of the solid and dotted fits, respectively. 
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Fig. 9. — The gravitational potential, in units of keV (see, e.g. 
Churazov et al. 2008), as a function of radius recovered from the 
X-ray analysis using Equation 4 (solid squares), with 1 a error bars 
from 100 Monte Carlo realisations of the method. For comparison 
wc show the gravitational potential recovered from the GC analysis 
of CR08 (stars), complete with 1 a error bars. The reference radius 
has been set (see text) such that the potential is zero at 14.2 kpc. 
The dashed line and associated 1 a confidence region shows the 
result of allowing for the abundance gradient in the analysis (see 
Section 4.1.1 for details). 
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+ C, (4) 



where C is an arbitrary constant and all other terms are 
defined as in Equation 3. This can be directly compared 
to the recovered potential from the dynamical mass anal- 
ysis, which has been determined in a Bayesian manner. 
To make the comparison, we need to assign a reference 
radius, at which the potential from both methods is set 
to zero. The choice of this radius is in fact arbitrary, but 
we have taken into account the properties of each profile 
in setting this radius to allow a useful comparison of the 
profiles. The shock regions in NGC 4636 extend to 100" 
(^7.8 kpc), and within this radius it is unlikely that the 
assumption of hydrostatic equilibrium is an adequate de- 
scription of the state of the gas. Therefore, we set the 
reference radius to be 14.2 kpc (~ 183"), which also takes 
into account the presence of close to 60 GCs within this 
radius, meaning that the potential from the dynamical 
mass analysis should be well-constrained. 

We compare the potential profiles recovered from each 
method in Figure 9. As expected from the comparisons 
of the total mass density shown in Figure 7, the poten- 
tial recovered from the X-ray analysis exceeds that from 
the GC analysis outside ~ 30 kpc. Again, the shape of 
the two profiles also differs outside a radius of 30 kpc; 
the X-ray gas and GCs trace the same underlying po- 
tential, so this is problematical. For comparison, we also 
show the results of allowing for the abundance gradient 
(shown as a dashed line), as explained in Section 4.1.1. 
Figure 9 shows that in the radial range ~100-300", the 
X-ray and dynamically derived potential profiles agree 
within the 1 a errors. Allowing for the metallicity gradi- 
ent in the X-ray data appears to improve the agreement 
at large radius, making the X-ray and dynamical profiles 
consistent within the quoted 1 a errors. However, care 
must be taken in the interpretation of these results, as 
the profiles have been normalised to equal zero at the 
same radius, and the agreement weakens if normalised 
elsewhere. The key point is that the shape of the X-ray 
and dynamically derived profiles agrees within the la er- 
rors over the radial range ~ 100" to ~ 300", but outside 
this radius, the gradient of the dynamical profile lessens 
with radius compared to the gradient of the X-ray po- 
tential profile. 

Churazov et al. (2008) explain in detail how directly 
comparing the potential recovered from each method 
lends some insight into the magnitude of any non-thermal 
pressure effects, as any non-thermal pressure support in 
the X-ray gas would lead to a smaller change in the X-ray 
potential compared to the dynamical potential. Hence, 
the gradient of a linear fit to Figure 10 would yield in- 
formation about the fractional contribution from non- 
thermal pressure support. However, we can immediately 
see that in this case there is not a simple linear expression 
linking the potentials, and the gradient increases with 
increasing potential. Following the prescription of Chu- 
razov et al. (2008), we cannot attribute the behaviour at 
large radius (potential) to the presence of non-thermal 
pressure support, as this would reduce the X-ray derived 
potential in relation to the dynamical potential. The ra- 
dial coverage of the GC data also limits the usefulness of 
this comparison for determining any non-thermal pres- 
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Fig. 10. — A direct comparison of the gravitational potential re- 
covered from the X-ray (y axis) and GC (x axis) analyses, evaluated 
at the radii of the dynamical profile (error bars). The grey confi- 
dence region and dashed lines show the effects of allowing for the 
abundance gradient in the analysis. The confidence region shows 
the 1 <r errors in the x direction, whereas the dashed lines show the 
1 a errors in the y direction. The solid line is for reference and 
shows a slope of 1. 

sure support component in the shocked region (< 100"). 



4.2.2. Mass-to-light ratio 

We can examine the central regions in detail, by com- 
paring the recovered if-band mass-to-light ratios from 
the two methods. The enclosed light profile was depro- 
jected by CR08 from a if-band surface brightness profile 
provided by Tom Jarrett (see CR08 for details) from the 
2MASS Large Galaxy Atlas (Jarrett et al. 2003). We use 
this enclosed light profile to determine the mass-to-light 
ratio from our X-ray mass profile. The comparison is 
shown in Figure 11, where the mass- to-light ratio from 
the X-ray analysis (solid line and associated confidence 
region) has been capped at the limit of the light profile 
data (23.5 kpc). Within 1 R e ff, the i-T-band mass-to- 
light ratio derived from the X-ray analysis decreases in- 
wards implying M*/Lk < 3; this is the region where the 
stellar mass component dominates over the dark matter 
component (Mamon & Lokas 2005). The mean stellar 
mass-to-light ratio in the if-band was found from 2dF- 
GRS and 2MASS data to be 0.73 M QiK /Lq iK assum- 
ing a Kennicutt IMF, and 1.32 Mq^k/Lq.k assuming a 
Salpeter IMF (Cole et al. 2001). This considered both 
early and late-type galaxies, but as the near-IR luminos- 
ity traces the older stellar population, it is reasonable 
to compare with this result. We also show the colour- 
dependent estimate of X-band stellar mass-to-light ratio 
from the prescription of Bell et al. (2003, see Section 
4.2 for more details). Humphrey et al. (2006) measured 
stellar mass-to-light ratios for 7 early-type galaxies, both 
from stellar population synthesis models and from mod- 
elling X-ray derived mass profiles with dark matter and 
stellar components. They find stellar mass-to-light ra- 
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Fig. 11. — The K-bsmd mass-to-light ratio from the GC analysis 
of CR08 (shown as open circles with error bars) and the ii"-band 
mass-to-light ratio derived from the X-ray analysis (shown as the 
solid line). The shaded region shows the combined errors from 
the light profile and the X-ray mass profile. The vertical dashed 
line shows the ii"-band R e f f of the galaxy from the 2MASS Large 
Galaxy Atlas (R e ff.K = 56.2", Jarrett et al. 2003), the horizontal 
dotted line shows the mean stellar mass-to-light ratio of Cole et al. 
(2001) (Kennicutt IMF) and the horizontal dashed line shows the 
colour-dependent M/L K of Bell ct al. (2003). 

tios ranging between ~0.5 and ~1.2 from the X-ray mass 
modelling. The stellar population models recover slightly 
higher values (~0.4 to ~1.9), depending on the assumed 
IMF. This indicates that the recovered i\~-band mass-to- 
light ratio in the central regions is consistent with pre- 
vious results, although the values at radii < 30" fall be- 
low the mean value of Cole et al. (2001) and the colour- 
dependent estimate of Bell et al. (2003). 

One possible explanation for the low stellar mass-to- 
light ratios seen in the central regions from the X-ray 
analysis is that in this region, we are underestimating the 
galaxy mass. We will explore this possibility in Section 
5. However, it is clear from Figure 11 that outside !R e ff 
there is a contributive mass component in addition to 
that expected from the stars alone, and this is predicted 
by both the X-ray and GC analyses. 

5. DISCUSSION 

We have shown in Figure 7 the derived total mass den- 
sity profile from our X-ray analysis, and the results of 
the dynamical analysis by CR08 of the GC system of 
NGC 4636. We have also compared the gravitational 
potential recovered from each method, and we note the 
following: 

1. Within ~10 kpc, the mass derived from the dynam- 
ical analysis of CR08 exceeds the mass recovered 
from the X-ray analysis. 

2. Between ~10 kpc and ~30 kpc, the profiles are 
consistent within the quoted la errors. 

3. Perhaps most crucially, outside ~30 kpc, the mass 
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recovered from the X-ray analysis significantly ex- 
ceeds that derived from the dynamical analysis. 

Considering the mass-to- light profile of the system sug- 
gests a significant dark matter component outside one ef- 
fective radius, independent of the analysis method. With 
the aim of understanding the behaviour in the inner and 
outer regions, we now review the observed discrepancies 
in terms of the key systematics of each analysis method. 

5.1. Anisotropy & CHASSIS 

The assumption that prevails within the current form 
of the dynamical analysis — the algorithm CHASSIS — 
is that of isotropy in phase space, and we should exam- 
ine the recovered dynamical density profile (Figure 7), 
in light of this assumption, or in particular, the assump- 
tion of velocity isotropy. At the outset, we note that 
deviation from isotropy in the true velocity space con- 
figuration of the system would urge CHASSIS to over- 
estimate the mass density, at radii where anisotropy pre- 
vails (Chakrabarty & Portegies Zwart 2004; Chakrabarty 
2006). If in reality, velocity anisotropy describes the 
phase space distribution from which the measured GC 
kinematic data are drawn, then the recovered mass den- 
sity distribution would be spuriously enhanced in am- 
plitude at these radii. In other words, the "true" mass 
density would be even lower than that indicated by the 
current dynamical estimates (see Figure 7). Therefore, 
invoking velocity anisotropy does not help to reconcile 
the X-ray and dynamical mass density profiles in the 
outer parts of NGC 4636. 

This however poses the question of whether the density 
distribution of CR08 is an overestimate, due to mistak- 
ing the velocity space configuration as isotropic, and if 
so, can we quantify how bad the assumption of isotropy 
is, given the measured kinematic data and our recovered 
density profile pi We need to find prob(a\{data}, p, K), 
where K is our state of background knowledge and a 
is a quantification of velocity anisotropy. For example, 
it could be parametrised in terms of the anisotropy pa- 
rameter (3. However, the maximum likelihood approach 
within CHASSIS calculates prob(f, p\ {data} , a = ao, K) , 
where / is the phase space density distribution that 
CHASSIS determines, along with p, and a is the value 
of a corresponding to isotropy in velocity space, within 
the adopted scheme of anisotropy parametrisation. In 
general, it is not possible to go from the calculated prob- 
ability to the required form. 

We can resort to an intermediate path by quantify- 
ing the probability of measuring a test statistic at least 
as extreme as the measured value of this statistic, given 
isotropy. This probability is referred to as a p-value. 
However, the p- value is a much maligned device, primar- 
ily because of the often neglected limitations of p- values 
and the subjectivity involved in establishing the accep- 
tance of a hypothesis. Also, the p-value is a probability 
defined on sample space, but it is more satisfying to work 
with an alternative obtained by considering the full pa- 
rameter space. 

We choose instead to employ the Bayesian evidence 
value or ev, details of which can be found in a well- 
written recent paper by Pereira et al. (2008) where a 
Fully Bayesian Significance Test (FBST) is advocated 
(Pereira & Stern 1999). Our null hypothesis H is that 



isotropy prevails in phase space. A brief synopsis of the 
FBST is presented in the Appendix, which in its full 
form requires the calculation of the Bayesian evidence 
value against Ho (ev), given by the integral of the pos- 
terior over the tangential set T. Here T comprises the 
mass density configurations p that correspond to poste- 
rior probability in excess of the posterior corresponding 
to p*, which in turn, is the point in p-space, that max- 
imises the posterior, while satisfying Hq, i.e. the max- 
imal p that stems from the assumption of isotropy. Fi- 
nally, ev is obtained as 1 — ev. The definition of FBST is 
that the test rejects Ho when ev is small. To ease our cal- 
culations, we view the integral over T in the conventional 
sense of treating probabilities, i.e. as the fractional num- 
ber of cases for which prob(p\{data}) > prob(p*). Here, 
the fraction is out of the total number N of recorded 
mass density distributions; since, one mass density dis- 
tribution is recorded for every iterative step, the fraction 
is calculated out of N, where N is the total number of 
iterative steps in a run of CHASSIS. 

Our simplification assumes that the N iterative steps 
cover the full parameter space. This may not be the 
case, though we need to remember that it is in propor- 
tion to the volume of the scanned parameter space that 
the volume of T is determined. In any case, the scanned 
range initiates with a seed (which has been established 
to be distant from the true configuration) and converges 
to the answer. We have also checked for the chain ex- 
tending to multiple times the burn-in period as well as 
it being well-mixed. Thus, we can bestow confidence on 
the recorded density distributions covering a substantial 
part of the parameter space. We find that for RUN I 
of CR08, ev~0.98, meaning our simpler version (over 
Pereira et. al's definition) of the Bayesian ev calcula- 
tion allows us to not reject H , i.e. not reject velocity 
isotropy as a valid assumption. In fact, following Stern 
(2000), we suggest that this high ev suggests "possibilis- 
tic support" in favour of the assumption of isotropy. 

5.2. The globular cluster system 

It is also worth noting the interesting features of the 
GC system of NGC 4636. A Chandra study of the 
level of association of GCs with LMXBs (Posson-Brown 
et al. 2009) has shown consistency with similar early-type 
galaxies (see Fabbiano 2006). The specific frequency of 
GCs in NGC 4636 has consistently been found to be high 
(~6-9, see discussion of Dirsch et al. 2005). Dirsch et al. 

(2005) also show that the radial distribution of GCs is 
shallower than the galaxy light within approximately 7'. 
The slope changes to be consistent with the galaxy light 
outside 7' for the red GCs; this occurs at ~9' for the 
blue population. Only 15 GCs are observed outside a 
projected radius of <~ 7.5' (Schuberth et al. 2006). 

Our statistical analysis of the effects of orbital 
anisotropy leads us to accept our null hypothesis of veloc- 
ity isotropy, indicating that the mass discrepancy at large 
radius is not the result of poorly handled anisotropy. The 
CHASSIS algorithm assumes the same distribution func- 
tion in phase space, and therefore assumes that each GC 
feels the same dark matter distribution. Schuberth et al. 

(2006) suggested that the break in the radial distribu- 
tion of GCs is inconsistent with NGC 4636 being em- 
bedded in a large dark matter halo. Figure 1 of CR08 
shows the distribution of GCs with measured velocities, 
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Fig. 12. — The three regions used for extracting spectra to com- 
pare the mass profile recovered in different directions from the 
galaxy centre. The regions are overlaid on the filtered ACIS-I 
events file, binned by 4x4. See Section 5.3.1 for details. 

and shows that outside 30kpc, the radial velocities are 
not symmetrical about zero. More extensive coverage of 
GC velocities in the radial range 30-60 kpc is required to 
thoroughly resolve this issue. 

5.3. X-ray systematics 

We can now examine the effects of the assumptions in- 
volved in the X-ray analysis, the most notable of which 
is the assumption of hydrostatic equilibrium. We also 
assess whether our choice of gas density model is appro- 
priate. 

5.3.1. Hydrostatic equilibrium 

The assumption of hydrostatic equilibrium is a pre- 
requisite in determining the mass profile from an X-ray 
analysis in the manner described by Equation 3. There 
is currently some controversy in this recent work 

by Diehl & Statler (2007) has proposed that the majority 
of early-type galaxy systems are not in hydrostatic equi- 
librium, which inevitably impacts the recovered masses. 
However, the results of work by Churazov et al. (2008) 
demonstrates that in mildly disturbed systems, agree- 
ment can be sought between X-ray and dynamically de- 
rived mass profiles, indicating that the assumption of 
hydrostatic equilibrium is indeed valid. By choosing to 
examine NGC 4636, we can assess the impact of any pos- 
sible departures from hydrostatic equilibrium in detail. 

To assess the assumption of hydrostatic equilibrium, 
we have compared the recovered mass profile in three 
'slices' through the coarse annuli in different directions 
(see Figure 12). We determine the mass profiles from 
these regions using the XSPEC PROJCT model, and it is 
envisaged that a disturbance in the gas affecting one of 
these slices will be visible in the recovered mass profile 
in comparison to the original analysis. This implictly 
assumes spherical symmetry, but will give an indication 
of the extent to which the mass profile is affected by 
looking at more localised regions, instead of averaging 
over a full annulus. 

Figure 13 shows the recovered deprojected tempera- 
ture profile, gas density profile and cumulative mass pro- 



file yielded from this analysis. The deprojection proce- 
dure in each case was more unstable than in our original 
procedure, due to the reduced number of counts in each 
spectrum, and we fixed the abundance to fit at a single 
value across all radii. The instability in the deprojection 
appears strongest in Region 3, where the penultimate 
temperature point is fitted very low, which if left in the 
fitting procedure, produces an unphysical decrease in the 
cumulative mass profile. We have therefore ignored this 
point in our smoothing spline fit. Figure 13 shows a 
broad consistency in the inner regions, although the re- 
sults from the 3 regions do differ, suggesting that as ex- 
pected, the central disturbances are affecting the mass 
profile under the assumption of hydrostatic equilibrium. 

However outside 300", the recovered profiles are all 
consistent with the original analysis, and do not agree 
with the lower mass from the dynamical estimate (shown 
as a grey confidence region) . This seems to indicate that 
the assumption of hydrostatic equilibrium is valid in the 
outskirts of this system. If localised disturbances in the 
gas were causing the assumption of hydrostatic equlib- 
rium to dramatically under or over-estimate the mass 
this would be visible in these profiles. This is a key re- 
sult; the local structure at large radius does lead to some 
small differences between the profiles, but the mass dis- 
crepancy between the X-ray and dynamical mass profiles 
is not the result of these structural differences. 

We can also examine the behaviour at small radius, 
where the X-ray mass is lower than the dynamically in- 
ferred mass. If we make the assumption that the dy- 
namically inferred mass is indeed the true mass, then we 
can postulate what the X-ray inferred mass is telling us 
about bulk motions in the gas. Ciotti & Pellegrini (2004) 
show that if an X-ray analysis is applied to a situation 
where the gas is not in hydrostatic equlibrium, the re- 
covered X-ray mass M est relates to the true mass M in 
the following way, 

M est = M + —v (5) 

where G is the gravitational constant, r is the radius 
and v describes the contribution from non-hydrostatic 
processes. This ignores the effects of pressure terms such 
as those from magnetic processes. Thus, in the inner re- 
gions of the profile where the dynamically inferred mass 
exceeds the X-ray mass profile, v must be negative; in 
the outer regions, v must be positive. In such a case, 
Ciotti & Pellegrini (2004) explain that in the central re- 
gion, gas must be inflowing, whereas the outer region 
must be outflowing. The next question is whether this 
situation can be physically maintained. It does appear 
so — Pellegrini & Ciotti (1998) show that these so-called 
•partial winds can exist, but it seems very unlikely that 
this system could host an outflow with enough velocity 
to affect the mass estimation. Mapping the gas proper- 
ties in detail outside > 30 kpc where we are reaching the 
group regime would help to settle this issue. 

Is it possible that the central region is hosting an in- 
flow? If this is the case, it must operate on small ra- 
dial scales, as the X-ray and dynamical mass profiles 
are consistent between ~ 150" and 400" indicating that 
hydrostatic equilibrium here is obeyed. However, this 
is unlikely to be a long-term inflow, due to the recent 
(3xl0 6 yr) outburst from the central AGN (Jones et al. 
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Fig. 13. — The results of the mass analysis using the three regions described in Figure 12. The leftmost column corresponds to the 
original analysis and shows the deprojected temperature profile (top), gas density profile (centre) and cumulative mass profile (bottom). 
The remaining three columns show the deprojected temperature profiles (top), gas density profiles (centre) and cumulative mass profiles 
(bottom) from regions 1, 2 and 3 (Figure 12) respectively. In each deprojected temperature profile, the solid line shows the smoothing 
spline fit to the profile. The solid lines in each gas density profile panel shown the /3-model fit to the data in that region. The dashed lines in 
the gas density profile panels show the same fit to the original data. The cumulative mass profiles are shown evaluated at the coarse radii. 
All errors shown are 1 a and come from 100 MC realisations (200 in the original analysis) of the procedure. The grey confidence region 
in the cumulative mass profile plots shows the cumulative mass profile recovered by the dynamical analysis of the GCs by Chakrabarty & 
Raychaudhury (2008). 



2002), and limits which can be placed on the total cool 
gas mass (Sage et al. 2007) suggest that the central re- 
gion is not cooling to form large deposits of cool gas. 
It is possible that the effect of the shocks was to push 
gas outwards, and it is now falling back. If the pres- 
ence of the shocks were affecting our spectral fits, they 
would raise the recovered temperature and density, and 
therefore would raise the mass, so the effect of excluding 
the shocks would be to lower the mass in the inner re- 
gions, thus worsening the discrepancy. We note that our 
X-ray solution which allows for the observed abundance 
gradient in fact reduces the recovered mass in the inner 
regions, due to the flattening of the gas density profile, 
hence worsening the agreement with the GC data in the 
inner regions (e.g. Figure 7). 

It is also possible that the mis-match in the inner re- 
gions is the result of additional non-thermal pressure sup- 
port in the gas, which could manifest from a variety of 
sources, such as rotation, the presence of magnetic fields, 
cosmic rays or a mixing of the gas with radio plasma. The 
presence of magnetic fields, cosmic rays or the mixing of 
the gas with radio plasma are all connected to the pres- 
ence of a central AGN. Although the galaxy hosts a cen- 



tral radio source, it only extends over the central 3 kpc 
(at 1.4 GHz; Jetha et al. 2007), but there is evidence of a 
previous, recent AGN outburst (Jones et al. 2002; Ohto 
et al. 2003) . Comparing the gravitational potential from 
the two methods can indicate the required non-thermal 
pressure contribution (Churazov et al. 2008), but the ra- 
dial coverage of the GC data prevents a comparison over 
the radius of interest. 

5.3.2. Gas density model 

Figure 3 shows that although a /3-model parameterisa- 
tion of the gas density profile performs well at ~100"- 
400", the central regions and the outer regions arc not 
well-described in this way. As the largest disagreement 
between the X-ray and dynamical profiles occurs at the 
largest radii, it is prudent to assess the implications of 
our model choice. We fit the gas density profile with a 
smoothing spline in linear-log space using the R PROJECT 
algorithm SMOOTH. SPLINE to better represent its shape 
and to allow for the local features in the profile to be 
incoporated into the analysis. The fit is shown in the 
top panel of Figure 14 (solid line), alongside the origi- 
nal /3-model (dashed line). The smoothing spline does 



14 



Johnson et al. 



co ° 
I o 

E 

CO 
CO 

D) ,- 
^§ 

o 




- 

_q_ CD 

CL O 

E 

zL o 

o 

CO 

'■4— ' CD 

C d 

CD i 
-i—' 

O CM 

Q_ v 



o 

CO 



Radius (arcsec) 

20 30 100 300 
— I 1 — | 



i — r 













— 1 1 1 1 1 1 


i 1 — 


• — H 



cgr°o 




30 100 300 

Radius (arcsec) 

Fig. 14. — Top panel: The fine binned gas density profile 
(shown as error bars) from Figure 3, with the associated /3-modcl fit 
(dashed line). The smoothing spline fit to the profile (in linear-log) 
space is shown as the solid line. Centre panel: The gravitational 
potential (in keV) recovered from the original analysis (shown as 
error bars, from Figure 9). The grey confidence region shows the 
dynamically derived potential and lo~ errors, again from Figure 
9. Open circles show the potential recovered from the smoothing 
spline fit to the gas density profile, shown in the top panel. Bot- 
tom panel: The cumulative mass profile from the original analysis 
(shown as error bars), reproduced from Figure 4. The grey confi- 
dence region shows the dynamically recovered mass profile of CR08 
(including la errors). The open circles show the mass profile re- 
sulting from the use of the smoothing spline fit to the gas density 
profile (shown in the top panel). 

not capture all the local features, but it performs well in 
representing the large radius behaviour. 

In Section 3.4, we noted the apparent flattening at 
large radius in the gas density profile reported by 
Trinchieri et al. (1994). The difficulty in making a con- 
clusive statement regarding the presence or not of a 
bump in our gas density profile lies in the inherent in- 
stability of the PROJCT model in such regions, where the 
surface brightness profile is at its flattest (see, for exam- 



ple Russell et al. 2008). This is a generic problem for de- 
projection schemes. In a physical system, M{< r) mono- 
tonically increases, but we can see from Figure 14 that 
allowing for a flattening in the gas density profile yields 
an unphysical mass profile. This is showing the limits of 
our deprojection, and to understand the gas properties 
at large radius requires mapping the gas properties to 
larger radius, beyond the scope of the current paper. 

6. CONCLUSIONS 

We present an X-ray mass analysis of the early-type 
galaxy NGC 4636 using Chandra data, under the as- 
sumptions of spherical symmetry and hydrostatic equi- 
librium. The integrity of the latter assumption has been 
questioned with reference to early-type galaxies (Diehl & 
Statler 2007), and it is because of the observed distur- 
bances in the gas in NGC 4636 (e.g. Jones et al. 2002; 
Ohto et al. 2003; O'Sullivan et al. 2005) that we chose 
to study this object, in an effort to assess the impact of 
this assumption on the recovered mass profile. We find 
that the treatment of the abundance gradient in the X- 
ray analysis can significantly affect the recovered mass 
profile at all radii. 

We have compared the X-ray mass density profile with 
that recovered from a dynamical analysis of the sys- 
tem's globular clusters (GCs), presented by Chakrabarty 
& Raychaudhury (2008). Inside 10 kpc, the dynamical 
mass estimate exceeds the X-ray mass estimate. The gas 
in this region is highly disturbed, and we postulate the 
cause of the disagreement to be a localised inflow of gas, 
or a contribution of non-thermal pressure support. 

The mass density profiles over the range ~10-30 kpc 
are consistent within 1 a, indicating that even in this 
highly disturbed system, the recovered X-ray mass is con- 
sonant with that recovered from an independent method 
over intermediate radii. However, outside 30 kpc, the X- 
ray mass estimate exceeds the dynamical mass estimate, 
by a factor of 4-5 times at its greatest disagreement. 
Examining the anisotropy of the GCs, we find no sta- 
tistical reason to reject our assumption of isotropy. The 
GC analysis is model- independent, so is not limited by 
the method, but the paucity of measured GC kinemat- 
ics outside 7.5' means that the success of this method at 
large radius is limited by the data. 

We test the assumption of hydrostatic equilibrium in 
our X-ray analysis, finding that local disturbances at 
large radii do not account for the observed discrepancy. 
At this radius, the group gas contribution is important 
in this system (O'Sullivan et al. 2005), and the overall 
state of the gas at this radius is uncertain. Mapping the 
X-ray properties to a larger radius using XMM-Newton 
would help to model the group emission, but is beyond 
the scope of the current paper. 

The X-ray and dynamical mass analysis methods both 
indicate the need for a dark matter halo in this system, 
and provide a useful comparison within 30 kpc. It is 
through the comparison of independent approaches that 
the most robust constraints will be placed on the mass 
distribution of early-type galaxies, but we conclude that 
the limiting factors in such a comparison to large radius 
(outside 30 kpc) are data quality in the case of the GC 
kinematics, knowledge of the overall state of the gas as we 
reach the group regime in the case of the X-ray analysis, 
or a combination of the two. 
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APPENDIX 

FULLY BAYESIAN SIGNIFICANCE TEST 

The Fully Bayesian Significance Test (Pereira & Stern 1999; Pereira et al. 2008) requires the computation of the 
Bayesian evidence value (ev), which we define here. Let prob(p) be the probability density function (pdf) over space g. 
Let the posterior probability of p given a measurement (represented by the data set {data}) is prob(p\{data}). Let p* 
be the point that maximises the posterior, while satisfying the null hypothesis H . Then, the evidence value against 
H is: 

ev = prob(p E T\{data}), (Al) 

where T is the tangential set, defined as: 

T = {pe g: prob(p\{data}) > prob(p*)}. (A2) 

Then ev can be written as: 

ev — J prob(p\{data})dp (A3) 
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